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An  Efficient,  Accurate  Numerical  Method  for  the 
Solution  of  a Poisson  Equation  on  a Sphere 


IMKODltTION 


The  need  for  efficient  and  accurate  methods  for  the  solution  of  Poisson-type 
equations  in  mathematical  physics  is  well  established.  Important  texts  in  numer- 
ical analysis  such  as  Matrix  Iterative  Analysis  (Varga^)  are  not  only  eloquent 
testimonies  to  the  magnitude  of  the  problem,  but  also  vivid  evidence  of  the  extent 
of  the  interest.  In  the  case  of  numerical  weather  prediction,  where  solutions  to 
the  Poisson  equation  are  required  in  daily  routine  operations,  it  is  paramount  that 
the  solution  procedure  be  efficient.  In  anticipation  of  the  development  of  global 
numerical  weather  prediction  models,  we  developed  earlier  an  efficient  shooting 
algorithm  for  the  solution  of  a discrete  Poisson  equation  on  the  surface  of  a sphere 
(Yee  ).  This  report  presents  the  most  recent  development  in  our  work  along  this 
line. 

We  seek  a numerical  solution  to  the  Poisson  equation  on  the  surface  of  a unit 
sphere. 


(Received  for  publication  28  October  1977) 

1.  Varga,  R.S.  (1962)  Matrix  Interative  Analysis,  Prentice-Hall,  Englewood 

Cliffs,  N.J.  ^ 

2.  Yee,  S.  Y.  K.  (1976)  An  efficient  method  for  a finite-difference  solution  of  the 

Poisson  equation  on  the  surface  of  a sphere,  J.  Comput.  Phys.  22:215-228 
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Here  0 < 0 < is  the  colatitude  and  0 < X«  2Jr  is  the  longitude.  Discretization  of 
this  equation  using  a five-point  centered-difference  operator  for  the  Laplacian 
gives  the  following  I X J linear  algebraic  system 


a.  u._^  -(aj  + b.)  u.  +b.  u.^^  -c.Ru.  =f. 


1 < i < I 


where 


u.  = u.  . 
I i.J 


u.  . = u(0.,  X.) 

1,3  ' r j'  ’ 


Ae  = ir/I 


AX  = 2:r/J  , 


sin  0.  , /- 
a.  = 1-1/2 

' sin  0. 
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Here  we  have  excluded  from  consideration  the  coordinate  singularities  at  the  poles 

3 

by  adopting  the  spherical  grid  system  proposed  by  Merilees.  We  have  also 
gridded  the  surface  of  the  sphere  with  I grid  points  along  a meridian  and  J grid 
points  along  a latitude.  Note  that  a solution  for  Eq.  (2)  exists  only  if  its  right- 
hand  side  satifies  the  compatibility  condition 


Y sine.  L.-O  . (3) 

i=l  j=l 

Even  then,  since  the  coefficient  matrix  of  Eq.  (2)  is  singular,  with  its  rank  one  less 
than  its  order  (a^  " *^I  “ solution  can  be  determined  only  to  within  an  add- 

itive constant.  This  is  consistent  with  the  fact  that,  for  the  lack  of  lateral  bound- 
aries on  a sphere,  a solution  to  Eq.  (1)  can  only  be  determined  to  within  an  additive 
constant. 

Equation  (2),  however,  is  prone  to  numerical  instability.  First,  the  coefficient 
matrix  of  Eq.  (2)  has  a spectral  radius  that  is  larger  than  unity.  Furthermore, 
because  of  the  convergence  of  meridians  on  a sphere,  the  condition  number  of  the 
system  increases  rapidly  with  increasing  spatial  resolution  of  the  computational 
grid.  These  instability  properties  severely  limit  the  usefulness  of  the  simple 
shooting  method  reported  earlier  to  relatively  coarse  grid  resolutions.  For  larger 
systems  resulting  from  finer  resolutions,  because  of  the  inherent  numerical  insta- 
bility mentioned  above,  direct  application  of  this  method  will  give  inaccurate  re- 
sults. It  is  the  purpose  of  this  report  to  demonstrate  that  this  numerical  instabil- 
ity problem  can  be  remedied  by  a two -pronged  approach.  The  instability  due  to  a 
spectral  radius  larger  than  unity  is  alleviated  by  the  use  of  a multiple  shooting 
technique,  while  that  due  to  a large  condition  number  is  overcome  by  the  use  of  a 
flexible  grid. 


2.  SYNOPSIS  OK  \IK  rilOl) 


To  view  the  shooting  method  described  here  from  a proper  perspective,  we 
shall  first  review  in  this  section  a classical  method  for  the  solution  of  differential 
equations.  It  will  then  be  easy  to  see  that  this  method  is  indeed  no  mere  than  the 
implementation,  in  the  discrete  domain,  of  a solution  method  in  differential  equa- 
tions. The  classical  approach  for  the  solution  of  a differential  equation  separates 


Merilees,  P.  E.  (1973)  Pseudospectral  approximation  applied  to  the  shallow 
water  equations.  Atmosphere  11:13-20. 


r 


I ' 


the  complete  solution  into  two  parts,  a particular  part  and  a homogeneous  part. 
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For  example  (Berg  and  McGregor  ),  for  the  boundary  value  problem 


^u  = f 


in  R 


(4) 


u = u on  B , 


it  is  often  easy  to  find  a particular  solution  Up  which  satisfies  only  the  differential 
equation 

Up  = f . 

Continuity  requirements  on  the  boundary  then  give  us  Up  values  of  Up  on  the  bound- 
ary. We  solve  next  for  Up  in  the  homogeneous  system 


V Up  = 0 in  R 


* * „ 

u,-  = u - u_  on  B . 

x1  P 


(5) 


The  complete  solution  for  Eq.  (4)  is  then  given  by 
u = Up  + Up  . 

Equation  (2)  is  a discrete  form  of  an  inhomogeneous  two-dimensional  linear  elliptic 
equation  applied  in  a closed  region  with  no  external  boundary.  It  may  be  looked 
upon  as  a discrete  boundary  value  problem  in  the  sense  that,  although  subject  to  no 
explicit  external  boundary  constraints,  it  must  satisfy  the  global  constraint  Eq.  (3). 
We,  therefore,  propose  to  solve  Eq.  (2)  in  two  steps.  First,  we  construct  a par- 
ticular solution  which,  except  at  certain  preselected  grid  latitudes  0^,  satisfies 
Eq.  (2).  One  way  of  constructing  such  a particular  solution  for  Eq.  (2)  is  to  assume 
arbitrary  values  for  the  components  of  Up  at  0p  We  may  then  compute  by  a march- 
ing procedure  a solution  of  the  equation  for  the  rest  of  the  region.  At  0.,  the  fore- 
ing  function  fj  computed  from  the  particular  solution  will  be  different,  in  general, 
from  the  given  fj.  The  particular  solution,  therefore,  does  not  satisfy  Eq.  (2)  at 
0j.  We  next  "correct"  these  discrepancies  by  adding  to  the  particular  solution  the 
homogeneous  solution  to  Eq.  (2).  Since  the  homogeneous  equation  has  only  J non- 
zero  components  (Tj  - Tj  ) in  its  forcing  function,  it  can  be  solved  with  less  com- 
puting effort  than  that  required  for  the  inhomogeneous  equation. 

4.  Berg,  P.  W. , and  McGregor,  J.  L,  (1966)  Elementary  Partial  Differential 
Equations,  Holden -Day,  San  Francisco,  CA. 
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shooting  method  computes  first  the  complete  solution  at  selected  grid  latitudes 

^ (the  "missing  initial  conditions").  The  complete  solution  for  the  remainder  of  the 

grid  is  then  obtained  by  a marching  procedure  (shooting).  The  key  idea  here  is 
i 

' that,  by  the  use  of  shooting,  we  avoid  most  of  the  transform  operations  needed  to 

compute  the  homogeneous  solution  for  the  entire  region. 

Although  the  idea  here  is  to  avoid  the  need  of  calculating  the  entire  homogen- 
, eous  solution,  looking  at  the  shooting  method  as  just  a way  of  implementing  a 

classical  method  in  differential  equation  is  useful  because  it  provides  us  with  the 
j insight  to  compute  easily  the  "missing  initial  conditions"  through  the  judicious 

' use  of  Fourier  transforms.  This  is  especially  desirable  in  the  case  of  multiple 

shooting  where  we  have  to  seek  the  "missing  initial  conditions"  at  a number  of 

I 

5 grid  latitudes.  We  emphasize  here  that  the  construction  of  a particular  solution 

I for  Eq.  (2)  is  accomplished  simply  by  marching  in  the  two-dimensional  physical 

. domain.  Details  of  computing  the  "missing  initial  conditions"  are  described  in 

Section  3 in  connection  with  the  multiple  shooting  algorithm. 

3.  TIIK  I SK  OF  MI  LTIPI.K  SHOOTING 

We  mentioned  in  Section  1 that  the  shooting  method  is  inaccurate  for  the 
numerical  solution  of  Eq.  (2)  when  I and  J are  large.  This  difficulty  can  be 
alleviated  somewhat  by  the  use  of  a multiple  shooting  technique,  a straightforward 
extension  of  the  approach  described  in  Section  2.  With  this  approach,  a sphere  is 
divided  into  subregions  by  pairs  of  internal  computational  boundaries.  A partic- 
ular solution  is  then  constructed  for  each  of  the  subregions.  Again  these  solutions 
satisfy  Eq.  (2)  except  at  certain  prespecified  6^.  They,  therefore,  constitute  a 
particular  solution  for  the  sphere  as  a whole.  To  this  particular  solution,  we 
now  add  the  homogeneous  solution  of  Eq.  (2)  to  form  the  complete  solution.  This 
time,  the  forcing  function  of  the  homogeneous  equation  has,  however,  a number 
of  pairs  of  nonzero  component  vectors  due  to  the  arbitrarily  picked  internal  com- 
putational boundary  conditions. 

For  pedagogic  purposes,  we  shall  detail  here  a method  for  obtaining  a homo- 
geneous solution  to  Eq.  (2).  The  method  makes  use  of  similarity  transforms  to 
reduce  first  the  block-tridiagonal  system  Eq.  (2)  to  a system  of  tridiagonal  sys- 
tems. The  homogeneous  solutions  to  these  reduced  systems  are  then  computed 

by  an  efficient  Gaussian  elimination  algorithm  (Varga^).  The  motivation  of  such 

2 

an  approach  has  been  discussed  earlier  (Yee  ). 

If  we  perform  a discrete  Fourier  transform  on  u.  and  T and  let 

w.  = P"^  u.  , g.  = 7.  (6) 

i 1 ’ 1 
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we  may  write  Eq.  (2)  in  a space-Fourier  domain  for  1 < i « I, 


ai  w._j  - (a.  + b.)  w.  + b.  w.^j  - c.  D w.  = g. 
where 


(7) 


i.  1 


=i.  1 


w. 

1 


i.  k 


=i,  k 


'i.K, 


=i.K 


,-l 


and  D = P RP,  K = J.  Here  the  diagonal  elements  of  the  diagonal  matrix  D are 
the  eigenvalues  dj^  of  R,  and  P is  a J X J orthogonal  matrix  whose  columns  are 
the  normalized  eigenvectors  associated  with  dj^  that  is, 


= 2(1  - cos  , 1 < k s K 


cos  , 

2 .1/2  ( (cos  j\)/2^/2 
k ■ V J / 

■ - j-k 

,1/2 


1/2' 


1 :<  k < (K/2)  - 1 , 

k = K/2  , 

(K/2)  + 1 < k < K - 1 
k = K . 


1 £ j £ J 


(8) 


Grouping  w.  by  components  and  defining 

/-i.A  Ai,k\ 


For  k < K,  is  nonsingular,  and  Eq.  (9)  can  be  solved  efficiently  by  Gaussian 
elinnination.  For  the  case  of  k = K,  since  d„  ir  " system 


is  singular  and  thus  has  no  unique  solution.  This  merely  confirms  that  the  finite- 
difference  system  Eq,  (2)  is  consistent  with  the  differential  system  Eq.  (1).  We 
may  pick  arbitrarily  a value  for  one  of  the  components  of  w.^  to  solve  Eq.  (10). 

Following  the  approach  outlined  in  Section  2,  we  construct  first,  by  marching, 
a solution  vector  w.'  which  satisfies 


Here  g^^  differs  from  gj^  only  in  as  many  components  as  there  are  internal  compu- 
tational boundaries  on  the  sphere.  For  example,  for  the  case  where  the  surface  of 
the  sphere  is  divided  into  two  halves  at  the  equator,  and  where  the  matching-in- 
the-middle  technique  has  been  used  to  halve  the  shooting  subrange  (Fox^),  we  have 


Fox,  L,  (1957)  The  Numerical  Solution  of  Two-Point  Boundary  Problems  in 
Ordinary  Differential  Equations,  Oxford  University  Press,  Fair  Lawn,  N.J 


; 1 


n. 

\ 

i 


li 

^ ' 


|i  '• 


I = i/i  . 


(12) 


Here  Aw!  , is  the  discrepancy  of  the  kth  component,  at  d.,  of  the  two  marching 
solutions  from  opposite  boundaries  of  a subregion.  Thus  if  is  a solution  to  the 
homogeneous  system  for  Eq.  (9), 


Tk\  = A4  . (13) 

then  the  complete  solution  to  Eq.  (9)  is 

• (14) 

Once  Wj^  is  available  for  all  k,  u^  may  then  be  computed  from  the  inverse  trans- 
form to  Eq.  (6), 


Ui  = Pw.  . (15) 

At  this  point,  it  is  appropriate  to  point  out  that  the  derivations  detailed  here 
serve  only  to  demonstrate  that  in  the  case  of  multiple  shooting,  the  "missing 
initial  conditions"  can  easily  be  obtained  by  solving  tridiagonal  systems  such  as 
Eq.  (13).  As  mentioned  earlier  in  Section  2,  for  reasons  of  computational 


12 


efficiency,  we  use  neither  Eq.  (12)  to  yield  nor  Eq.  (15)  to  obtain  iT.  In 

practice,  may  be  computed  by  the  use  of  a much  more  efficient  shooting  algo- 
rithm (Yee^). 

Test  computations,  using  the  multiple  shooting  technique,  have  been  made  for 
system  (2)  for  various  n,  the  number  of  shooting  subregions,  and  for  several 

values  of  I and  J,  The  design  of  the  numerical  tests  is  the  same  as  that  reported 
2 

in  Yee.  Namely,  we  create  a set  of  true  solution  v^  ^ for  Eq.  (2)  by  computing 
f.  . from 


1,1 


= a.  V.  , . - (a.  + b.  + 2c.)  v.  . + b,  v. . , . + c.(v.  . , + v.  . , , ) 
1 1-1.3  II  1 1,3  I ‘+1,3  ‘ ‘,3-1  ‘,3+1 


Here  values  for  v^  ^ are  obtained  from  a random  number  generator  and  have  been 
subjected  to  the  constraint  E C v.  ^ = 0 so  that  f.  ^ satisfies  the  compatibility  con- 
dition Eq.  (3).  With  the  forcing  function  T ^ computed  in  this  manner,  a normal- 
ized error  norm  defined  by 


IIe|L  = IIu  - v|!  /||v|| 


may  then  be  considered  as  a measure  of  the  accuracy  of  the  numerical  procedure. 
The  number  of  digits  of  accuracy  in  u is  then  given  by 


Z = -logjo  IIEII2  . 

Table  1 gives  sample  results  for  a 64  X 64  grid  in  terms  of  Z,  ||  EII2  and 
I E(Max)l , the  magnitude  of  the  maximum  error  over  the  entire  computational 
region.  Here  n is  the  number  of  subregions  into  which  tha  surface  of  the  sphere 
is  divided  and  m is  the  number  of  grid  latitudes  within  each  subregion.  Thus  for 
n = 1,  m = 32,  for  example,  the  entire  surface  of  the  sphere  is  t.reated  as  one 
region,  and  shooting  is  conducted  from  both  poles  over  a shooting  range  of  32  grid 
latitudes  toward  the  equate”.  Inspection  of  Table  1 gives  us  the  impression  that 
for  the  case  of  a 64  X 64  s„  n this  method  is  useful  only  when  the  sphere  is 
divided  into  a relatively  lar^e  ,mber  of  latitude  bands.  For  example,  for  the 
case  of  n = 4,  we  have  Z = 1.6,  that  is,  less  than  two-digit  accuracy  in  our  results. 
Careful  study  of  the  distribution  of  the  error  fields  suggests,  however,  that  due  to 
the  convergence  of  meridian  on  a sphere,  the  largest  local  errors  invariably 
occur  in  the  computational  subregions  containing  the  polar  caps.  Furthermore, 
these  errors  are  considerably  larger  in  magnitude  than  those  in  other  subregions. 
This  phenomenon  is  well  borne  out  by  Figure  1 which  depicts  typical  error  dis- 
tributions as  a function  of  6.  for  the  cases  of  n = 4 and  n = 8.  We  see  from  this 
figure  that,  except  for  the  subregions  containing  the  polar  caps,  we  have  for  the 
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Table  1.  Error  Norma  in  the  Computed  Solution*  as  a Function  of  m 


m 

n 

Z 

IIEII3 

|E(Max)I 

2 

16 

9.  6 

2.  2 X 10"^° 

1.  0 X 10'® 

4 

8 

6.4 

3.  6 X 10"’^ 

2.  8 X 10'® 

8 

1.  6 

2.  3 X 10'^ 

2.  0 X 10'^ 

16 

-6.  2 

1.7X10® 

2.  4 X 10'^ 

32 

1 

-21.  4 

2 1 

2,  2 X 10 

3. 4 X 10^^ 

Z Number  of  digits  of  accuracy  in  computed  solution  u. 
* For  a 64  X 64  grid  with  constant 


GRID  LATITUDE 

Figure  1.  Effect  of  Subdivision,  Standard  Grid 


case  of  n ^ 4 more  than  two-digit  accuracy  in  our  results.  Comparison  of  the 
error  distributions  of  these  two  cases  also  suggests  that  increasing  the  number  of 
subregions  is  not  an  effective  way  to  reduce  the  deleterious  effect  due  to  the  con- 
vergence of  meridians.  Other  more  effective  means  must  be  sought  if  we  want  to 
use  a grid  with  high  resolution. 


1.  TIIK  I SK  OK  FLKMlll.K  OHIO 

It  is  easy  to  see  that  the  numerical  instability  inherent  in  the  solution  of 
Eq.  (2)  is  due  in  part  tc  the  convergence  of  meridians  on  a sphere.  From  the 
definition  of  Cj^  in  Eq.  (2), 


c.  = 

i 


^2  .2 
sin  0. 


2 

we  see  that  c,  ->■  oo  as  sin  6.  0.  Thus  if  is  not  a function  of  0..  increased 

1 i 1’ 

resolution  in  A0  tends  to  destabilize  the  solution  through  the  rapid  increase  of  the 
condition  number  of  the  coefficient  matrix  of  Eq.  (2).  We  shall  demonstrate  in  this 
section  how  this  difficulty  can  be  circumvented  by  the  use  of  a flexible  finite- 
difference  scheme. 

Consider  the  following  finite -difference  approximation  to  the  second  term  of 
Eq.  (1), 


1 


sin  0. 


0(13?) 


(16) 


where  3-  = From  the  numerical  stability  point  of  view,  it  is  desirable  to  keep 

^ ^ 2 

the  value  of  (3-  sin  0.)  from  approaching  zero.  The  simplest  way  of  achieving 
this  is  of  course  to  set  3j  = 1/sin  0..  This  is,  however,  impractical  because  the 
points  (i,  j-3j),  (i,  j+3j)  in  Eq.  (16)  will  not  in  general  coincide  v.'ith  the  regular 
grid  points.  To  insure  that  the  points  (i,  j±3j)  coincide  with  the  grid  points,  we 
must  have  3j  take  on  only  multiple  values  of  AA,  that  is. 


3i  = hj  A^ 

where  h.  are  positive  integers.  Since  considerations  of  maintaining  a uniform 
increment  in  actual  distance  on  a sphere  give  us 


15 


h'  = i — ^ — 

i I sm  0. 


we  therefore  simply  pick  as  the  integer  values  which  are  closest  to  h^.  Linder 
these  restrictions,  we  may  now  write  Eq.  (IG)  as 


1 3^u 

sin^  0 


'1.3 


- . 2 .^2  . 2 ^ 
h.  AA.  sin  6. 
i 1 


(U; 


3 -hi 


2u.  . + u.  ...  ) + — 
1,3  i,  3+hi  .....2 


in  n. 

I 


C)(h^  AA“) 


(17) 


Table  2 gives  the  values  of  h.  determined  in  this  manner  for  the  case  of  a 
64  X G4  grid.  Thus,  for  example,  for  i = 1,  we  have  h.  = h^,^  = 10,  indicating  that 
at  the  grid  latitudes  next  to  the  poles,  our  spatial  increments  in  A are  taken  over 
a distance  of  16  aA.  In  cases  where  I is  very  large,  precaution  should  be  taken  to 
insure  that  the  values  of  h^  near  the  poles  are  small  compared  to  J. 


Table  2.  Spatial  Increments  in  A as  Func- 
tion of  Latitude  for  a 64  X 64  Grid 


1 

h.,  h„_  . 
i'  65-1 

(unit;  aA) 

1 

16 

2 

14 

3 

CO 

4 

6 

5 

5 

6 

4 

CO 

3 

9,  . . . , 15 

to 

16,  . . .,  32 

1 

With  the  use  of  this  differencing  scheme,  Eq.  (2)  becomes 


aj  U..J  - (aj  + b.)  u.  + b.  u.^j 


c. 

R . u . = T. 

,2  i i i 


(18) 


i 

i 

i 

•1 

I 

1 

\ 

f 

1 


where  are  band-tridiagonal  J X J matrices  with  their  nonzero  off-diagonal  ele- 
ments separated  from  their  diagonal  elements  by  h.  zeros.  The  element  r.  . of 
R.  are  given  by 

r,  . = 2 , 1 = j , 

t,  J ' ' 

t = j ± h. 

= 0 , otherwise. 

Note  that  when  h.  = 1,  R.  become  the  matrix  R in  Eq.  (2).  Furthermore,  just  as 
R is  orthogonally  similar  to  a real  diagonal  matrix  D,  R.  are  orthogonally  similar 
to  real  diagonal  matrices  Dj, 

D.  = P'^  R.  P . 

Here  the  elements  of  P are  given  in  Eq.  (8)  and  the  diagonal  elements  of  are 
the  eigenvalues  d.  , of  R.,  that  is,  for  a given  h., 

If  K 1 1 

d^  it  = 2(1  - cos  hj  X^)  , 1 < k < K . 

Thus  with  the  exception  that  R^  are  now  latitude  dependent  and  we  must  compute 
the  eigenvalues  d.  for  different  h.,  Eq.  (18)  may  be  solved  in  exactly  the  same 
manner  Eq.  (2)  is  solved. 

Table  3 tabulates  sample  results  when  Eq.  (18)  instead  of  Eq.  (2)  is  used  as 
the  finite -difference  analog  of  Eq.  (1).  Again  the  results  are  for  a 64  X 64  grid 
and  the  accuracy  of  the  results  is  measured  in  terms  of  error  norms.  For  exam- 
ple, for  the  case  of  n = 4,  we  now  have  Z = 9.9,  almost  10-digit  accuracy  in  our 
numerical  results.  Comparison  of  Tables  1 and  3 clearly  shows  that  Eq.  (18)  is 
numerically  far  more  stable  than  Eq.  (2).  Note  that  for  a given  grid  resolution, 
we  may  choose  from  tables  such  as  these  values  of  n to  suit  our  needs,  depending 
on  the  desired  accuracy  in  our  results  and  the  number  of  digits  of  accuracy  carried 
by  the  computer.  (Our  CDC  6600  carries  roughly  15  digits.) 

A typical  error  distribution  as  a function  of  latitude  for  the  results  of  Eq.  (18) 
is  given  in  Figure  2 for  a case  of  n = 4.  Also  plotted  in  Figure  2 is  the  error  curve 
for  the  case  of  n = 4 from  Figure  1.  Comparison  of  these  two  curves  clearly 
demonstrates  that  the  difficulty  we  had  with  the  use  of  Eq.  (2)  has  now  disappeared. 


1 < j < J : 


. 1 =5  j - h. 

, h.  < j < (J  - h.) 

, (J  - h.  + 1)  < j « J 


I 
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Table  3.  Error  Norms  In  the  Computed  Solution*  as  a Function  of  m 


m 

n 

z 

IIEII3 

1 E(Max)| 

2 

11.0 

9.  4 X 10'^^ 

1. 1 X id'^^ 

4 

11.0 

9.  8 X 10'^^ 

2. 6 X 10'^^ 

8 

9.  9 

1.3  X 10”^° 

1.  5 X 10"^ 

16 

6.  5 

3.  3 X 10"’^ 

5.0  X 10‘® 

32 

1^1 

2.  2 

6.3X10"^ 

7.0  X 10*^ 

Z Number  of  digits  of  accuracy  in  computed  solution  u. 

* For  a 64  X 64  grid  with  flexible  finite-differencing  in  X. 
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3.  I OM  I.I  IMN(.  KKMVKKS 

The  following  observations  are  given  as  concluding  remarks  for  the  work 
' reported  here: 

‘ (1)  The  Hexible  finite -difference  scheme  described  in  Section  4 is  valuable 

- not  only  as  a means  of  stabilizing  the  finite-difference  shooting  solution  for  a 

Poisson  equation  on  a sphere,  but  it  is  valuable  also  for  the  finite -difference  for 
formulation  of  initial  value  problems  on  a sphere,  because  it  enables  us  to  relax  the 
severe  restriction  on  the  time  increments  near  the  poles  in  our  time  integrations. 

^ The  savings  in  computer  time  resulting  from  this  can  be  substantial. 

' (2)  Other  than  asserting  that  it  is  efficient,  we  have  not  discussed  here  the 

computational  efficiency  of  this  method  in  relation  to  other  direct  methods.  That 

; this  method  is  more  efficient  than  many  other  known  methods  has  been  demon- 

J 2 

' ^ strated  in  an  earlier  report  (Yee  ). 

: ^ (3)  As  a finite -difference  analog  to  Eq.  (1),  the  truncation  error  term  is 

f • larger  in  Eq.  (18)  than  in  Eq.  (2)  by  a factor  of  1/sin^  0..  Equation  (18)  is  thus 

, less  desirable  than  Eq.  (2)  from  this  point  of  view.  This  argument  in  favor  of 

; Eq.  (2)  is  lessened  somewhat  if  we  recall  that  other  sources  of  error  in  u (obser- 

i vational,  round-off)  are  all  inevitably  influenced  by  the  same  factor  1/sin^  0.. 

I (4)  Since  the  net  effect  of  the  flexible  grid  is  to  remove  the  effect  of  the  con- 

vergence of  meridians  on  a sphere,  the  results  in  Section  4 indirectly  indicate  that 
J the  multiple  shooting  technique,  without  incorporating  the  flexible  grid,  is  a viable 

tool  for  solving  the  Poisson  equation  in  geometries  such  as  rectangles  and  channels, 
j (5)  Since  Fourier  transforms  are  now  used  only  at  a few  selected  grid  lati- 

[ tudes,  we  may,  if  the  situation  warrants,  elect  not  to  use  fast  Fourier  transform 

for  increased  efficiency.  This  will  remove  the  severe  restriction  on  J.  J = 2*^, 
j where  k is  a positive  Integer. 
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